 
      module exbdrift
!---------------------------------------------------------------------- 
! description: calculates ExB drift velocities UE,VE,WE
!     uses the electric field which is calculated in module efield
!     on a regular magnetic grid (MLT deg/ mag. latitude)
!
!     0. initilize called before time-loop (exbdrift_init)
!     every timestep and every processor
!     1. map from magn. grid to geographic grid WACCM (map_mag2geo)
!     2. rotate e-field (rot_efield)
!     3. calculate ExB drift velocities  ue,ve,we [m/s] (iondrift)
!
! input ed1,ed2,    &    ! zonal/meridional elect. field [V/m]
!      nmlon,nmlat, &    ! dimension of mag. grid 
!      dlatm,dlonm, &    ! grid spacing of mag. grid 
!      ylatm,ylonm       ! magnetic MLT deg./latitudes (deg)
!
! ExB electromagnetic drift velocity [m/s] (east,north,upward)
! These are output to the physics buffer.
!      real(r8) ::  ue,ve,we
!
! notes: 
! - assume regular magnetic grid for the e-field (0:360 MLT/0:90) -> mapping
!
! Author: A. Maute Dec 2003
!         B. Foster adding physics buffer Feb, 2004.
!---------------------------------------------------------------------- 

      use shr_kind_mod,  only: r8 => shr_kind_r8
      use physconst,      only: pi
      use phys_buffer ,  only: pbuf_size_max, pbuf_fld, pbuf_add
      use ppgrid      ,  only: pcols      ! np. of atmospheric columns
      use cam_logfile,   only: iulog

      use efield, only : &  ! inputs from efield module
        ed1, ed2,     &     ! global zonal/meridional elect. field [V/m]
        potent,       &     ! electric potential [V]
        nmlon, nmlat, &     ! dimension of mag. grid 
        dlatm, dlonm, &     ! grid spacing of mag. grid 
        ylatm, ylonm        ! magnetic longitudes,latitudes (deg) (0:nmlat),(0:nmlon)

      implicit none

      private

      save

!---------------------------------------------------------------------- 
! Public interfaces:
!---------------------------------------------------------------------- 
      public :: exbdrift_init
      public :: exbdrift_register       ! register drift velocities with pbuf 
      public :: get_exbdrift

!---------------------------------------------------------------------- 
! Indices to drift velocities in physics buffer:
!---------------------------------------------------------------------- 
      integer  :: ndx_ue, ndx_ve, ndx_we
      real(r8) :: rtd, hr2d
      logical, parameter :: debug =.true.

      contains

      subroutine exbdrift_init  
!-----------------------------------------------------------------------
! Purpose: Prepare fields for histories.
!
! Method: 
!
! Author: A. Maute Dec 2003  am 12/30/03    
!-----------------------------------------------------------------------

   use cam_history,  only: addfld, phys_decomp, add_default ! for history saves

!-----------------------------------------------------------------------
! Add mag field output to master field list:
!-----------------------------------------------------------------------
   call addfld('UE      ','m/s     ',1,'I', 'Zonal ExB drift',phys_decomp)
   call addfld('VE      ','m/s     ',1,'I', 'Meridional ExB drift',phys_decomp)
!  call addfld('WE      ','m/s     ',1,'I', 'Vertical ExB drift',phys_decomp)
!  call addfld('MLT     ','hrs     ',1,'I', 'mag.local time',phys_decomp)
!  call addfld('EF_EAST ','V/m     ',1,'I', 'eastward electric field',phys_decomp)
!  call addfld('EF_WEST ','V/m     ',1,'I', 'northward electric field',phys_decomp)
!  call addfld('EF_UP   ','V/m     ',1,'I', 'upward electric field',phys_decomp)
!  call addfld('EF1_MAP ','V/m     ',1,'I', 'map. mag. eastward ef',phys_decomp)
!  call addfld('EF2_MAP ','V/m     ',1,'I', 'map. mag. northward ef',phys_decomp)
   call addfld('EPOTEN  ','V       ',1,'I', 'Electric Potential',phys_decomp)

!-----------------------------------------------------------------------
!  Write these fields to history by default:
!-----------------------------------------------------------------------
   call add_default ('UE' , 1, ' ')
   call add_default ('VE' , 1, ' ')	
!  call add_default ('WE' , 1, ' ')	
!  call add_default ('MLT', 1, ' ')
!  call add_default ('EF_EAST' , 1, ' ')
!  call add_default ('EF_WEST' , 1, ' ')	
!  call add_default ('EF_UP'   , 1, ' ')
!  call add_default ('EF1_MAP' , 1, ' ')	 
!  call add_default ('EF2_MAP' , 1, ' ')
   call add_default ('EPOTEN ' , 1, ' ')

      rtd  = 180._r8/pi     ! radians to degrees
      hr2d = 360._r8/24._r8

      end subroutine exbdrift_init

      subroutine exbdrift_register
!-----------------------------------------------------------------------
! Register drift velocity outputs with physics buffer:
!
! Ion velocities are 2d fields in pbuf (no vertical dimension),
!   so fdim = mdim = ldim = 1.
! Indices are saved as module data above.
!-----------------------------------------------------------------------

      call pbuf_add('UE', 'physpkg', 1, 1, 1, ndx_ue ) ! zonal
      call pbuf_add('VE', 'physpkg', 1, 1, 1, ndx_ve ) ! meridional
      call pbuf_add('WE', 'physpkg', 1, 1, 1, ndx_we ) ! vertical

      end subroutine exbdrift_register

      subroutine get_exbdrift( lchnk, ncol, pbuf )
!-----------------------------------------------------------------------
! Purpose: calculate ion drift velocities  [m/s]
!
! Method: v = E x B/B^2
!
! Author: A. Maute Dec 2003  am 12/30/03    
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ncol    ! np. of atmospheric columns
      integer, intent(in) :: lchnk   ! chunk identifier
      type(pbuf_fld), intent(inout) :: pbuf(pbuf_size_max) ! physics buffer

!-----------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------
      real(r8) :: ed1_geo(ncol), & ! electric field on geographic grid [V/m]
                  ed2_geo(ncol), &
                  epot_geo(ncol)   ! electric potential on geographic grid
      real(r8) :: elfld(3,ncol)    ! electric field in geog. direction on geographic grid [V/m]
      real(r8) :: mlt(ncol)        ! mag.local time of WACCM geo. grid point

!-----------------------------------------------------------------------
! calculate the magn. local time of WACCM geogr. grid points
!-----------------------------------------------------------------------
      call cal_mlt( mlt, lchnk, ncol )
!-----------------------------------------------------------------------
! map the electricfield from regular mag.grid to WACCM grid 
!-----------------------------------------------------------------------
      call map_mag2geo( mlt, lchnk, ncol, ed1_geo, ed2_geo, epot_geo )
!-----------------------------------------------------------------------
! rotate the electric field in geographic direction 
!-----------------------------------------------------------------------
      call rot_efield( lchnk, ncol, ed1_geo, ed2_geo, elfld )
!-----------------------------------------------------------------------
! calculate ExB drift velocities  ue,ve,we [m/s]
!-----------------------------------------------------------------------
      call iondrift( lchnk, ncol, elfld, pbuf )

      end subroutine get_exbdrift

      subroutine map_mag2geo( mlt, lchnk, pcol, ed1_geo, ed2_geo, epot_geo )
!-----------------------------------------------------------------------
! Purpose: map electric field from regular magnetic grid 
!    to WACCM physics geographic grid
!
! Method:  bilinear interpolation
!    assumptions: magnetic grid regular (MLT[deg]=0:360 deg/lat=0:180 deg)
!    ed1_geo(:,:) = Ed1_geo <- mapping of - 1/[R cos lam_m] d PHI/d phi_m
!    ed2_geo(:,:) = Ed2_geo <- mapping of 1/R d PHI/d lam_m/ sinIm
!
! Author: A.Maute Dec 2003 am 12/17/03 
!-----------------------------------------------------------------------

   use phys_grid,   only:  get_lat_p, get_lon_p ! latitudes/longitudes of current column (deg)
   use mo_apex,     only:  alatm                ! apex mag latitude at each geographic grid point (radians)
   use cam_history, only:  outfld

!-----------------------------------------------------------------------
! dummy arguments
!-----------------------------------------------------------------------
   integer, intent(in)   :: pcol      ! np. of atmospheric columns
   integer, intent(in)   :: lchnk     ! chunk identifier
   real(r8), intent(in)  :: mlt(pcol) ! mag.local time of WACCM geo. grid point
   real(r8), intent(out) :: &
     ed1_geo(pcol), &                 ! electric field on geog. grid
     ed2_geo(pcol), &                 ! electric field on geog. grid
     epot_geo(pcol)                   ! electric potential on geog grid

!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
     integer  :: i, iphi1, iphi2, ilam1, ilam2
     real(r8) :: t, u, collat, collon


     do i = 1,pcol
       collat = alatm(i,lchnk)*rtd      ! mag lats (deg) [-90:90]
       collon = mlt(i)*hr2d             ! mlt (deg) [0:360]
       iphi1 = int(collon/dlonm)	! indices of the 4 surrounding points (regular mag. grid
       iphi2 = iphi1 + 1	        ! with dlonm/dlatm deg grid spacing)
       ilam1 = int( (collat + 90._r8)/dlatm )
       ilam2 = ilam1 + 1
       if(iphi1 == nmlon ) then  ! boundaries 
           iphi1 = nmlon -1
           iphi2 = nmlon
       end if
       if(ilam1 == nmlat ) then 
           ilam1 = nmlat-1
           ilam2 = nmlat
       end if
       if(collon == 360._r8) then
          collon= 0._r8
          iphi1 = 0
          iphi2 = 1
       end if

       t = (collon - ylonm(iphi1))/(ylonm(iphi2) - ylonm(iphi1))
       u = (collat + 90._r8 - ylatm(ilam1))/(ylatm(ilam2) - ylatm(ilam1))
       ed1_geo(i) = (1._r8 - t)*(1._r8 - u)*ed1(iphi1,ilam1) + &
        	    t*(1._r8 - u)*	  ed1(iphi2,ilam1) + &
        	    t*u*	  ed1(iphi2,ilam2) + &
        	    (1._r8 - t)*u*	  ed1(iphi1,ilam2) 
       ed2_geo(i) = (1._r8 - t)*(1._r8 - u)*ed2(iphi1,ilam1) + &
        	    t*(1._r8 - u)*	  ed2(iphi2,ilam1) + &
        	    t*u*	  ed2(iphi2,ilam2) + &
        	    (1._r8 - t)*u*	  ed2(iphi1,ilam2) 
       epot_geo(i)= (1._r8 - t)*(1._r8 - u)*potent(iphi1,ilam1) + &
        	    t*(1._r8 - u)*	  potent(iphi2,ilam1) + &
        	    t*u*	  potent(iphi2,ilam2) + &
        	    (1._r8 - t)*u*	  potent(iphi1,ilam2) 
     end do ! i = 1,pcol

!     call outfld( 'EF1_MAP', ed1_geo, pcol, lchnk)
!     call outfld( 'EF2_MAP', ed2_geo, pcol, lchnk)
      call outfld( 'EPOTEN', epot_geo, pcol, lchnk)

     end subroutine map_mag2geo

   subroutine rot_efield( lchnk, pcol, ed1_geo, ed2_geo, elfld )
!-----------------------------------------------------------------------
! Purpose: rotate the electric field to get geographic east and westward direction
!
! Method: Richmond: J Geomag. Geoelectr. [1995] eqn. (4.5)
!      rotation of the electric field to get geog. eastward and downward/equatorward
!      E(k) = d_1(k)*ed1_geo + d_2(k)*ed2_geo  for k = 1,3
!
! Author: A. Maute Dec 2003  am 12/17/03 
!-----------------------------------------------------------------------

   use mo_apex,     only: d1vec, d2vec ! base vectors (3,pcols,begchunk:endchunk)
   use cam_history, only: outfld
      
!-----------------------------------------------------------------------
! dummy arguments
!-----------------------------------------------------------------------
   integer, intent(in)   :: pcol                         ! np. of atmospheric columns
   integer, intent(in)   :: lchnk                        ! chunk identifier
   real(r8), intent(in ) :: ed1_geo(pcol),ed2_geo(pcol)  ! electric field on geog. grid
   real(r8), intent(out) :: elfld(3,pcol)                ! electric field on geog. grid geog. direction

!-----------------------------------------------------------------------
! local
!-----------------------------------------------------------------------
   integer :: i, k

   do i = 1,pcol
     do k = 1,3
       elfld(k,i) = ed1_geo(i)*d1vec(k,i,lchnk) + & 
                    ed2_geo(i)*d2vec(k,i,lchnk)
     end do
   end do

!   call outfld( 'EF_EAST', elfld(1,:), pcol, lchnk)
!   call outfld( 'EF_WEST', elfld(2,:), pcol, lchnk)
!   call outfld( 'EF_UP', elfld(3,:), pcol, lchnk)

   end subroutine rot_efield

   subroutine iondrift( lchnk, pcol, elfld, pbuf )
!-----------------------------------------------------------------------
! Purpose: calculate ion drift velocity
!
! Method: v_i = ExB/B^2
!      for high altitudes where collisionfrequency mue_in << omega_i and
!      mue_en,vert << omega_e
!      v_i,vert = v_e,vert = v_E
!      B magnetic field component from apex code
!         bnorth northward gauss
!         beast  eastward  gauss
!         bdown  downward  gauss -> for upward component -bdown
!         bmag   magnitude gauss
!
! Author: A. Maute Dec 2003  am 12/17/03 
!-----------------------------------------------------------------------

   use mo_apex,     only: beast, bnorth, bdown, bmag ! component of B-field, |B| [Gauss]
   use cam_history, only: outfld
   use phys_grid,   only: get_lat_p, get_lon_p       ! latitudes/longitudes of current column (deg)

!-----------------------------------------------------------------------
! dummy arguments
!-----------------------------------------------------------------------
   integer, intent(in)  :: pcol           ! np. of atmospheric columns
   integer, intent(in)  :: lchnk          ! chunk id
   real(r8), intent(in) :: elfld(3,pcol)  ! electric field on geog. grid geog. direction
!-----------------------------------------------------------------------
! Ion velocities are saved in physics buffer:
!-----------------------------------------------------------------------
   type(pbuf_fld), intent(inout) :: pbuf(pbuf_size_max) ! physics buffer

!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
   integer  :: i
   real(r8) :: fac
   real(r8), pointer :: ue(:), ve(:), we(:)
 
!-----------------------------------------------------------------------
! Set local pointers to respective fields in pbuf:
!-----------------------------------------------------------------------
   ue => pbuf(ndx_ue)%fld_ptr(1,1:pcols,1,lchnk,1)
   ve => pbuf(ndx_ve)%fld_ptr(1,1:pcols,1,lchnk,1)
   we => pbuf(ndx_we)%fld_ptr(1,1:pcols,1,lchnk,1)

   do i = 1,pcol    ! number of columns in each chunk
!-----------------------------------------------------------------------
! Magnitude of magnetic field bmag [nT] is from apex module. 
!-----------------------------------------------------------------------
      fac = 1.e9_r8/bmag(i,lchnk)**2           ! nT to T (nT in denominator)
      ue(i) =  -(elfld(2,i)*bdown(i,lchnk) + elfld(3,i)*bnorth(i,lchnk))
      ve(i) =   elfld(3,i)*beast(i,lchnk) + elfld(1,i)*bdown(i,lchnk)
!     we(i) =   elfld(1,i)*bnorth(i,lchnk) - elfld(2,i)*beast(i,lchnk)
      ue(i) = ue(i)*fac
      ve(i) = ve(i)*fac
!     we(i) = we(i)*fac
   end do ! i = 1,pcol	

   call outfld( 'UE', ue, pcols, lchnk)
   call outfld( 'VE', ve, pcols, lchnk)
!  call outfld( 'BNORTH', bnorth(:pcol,lchnk)*1.e-5_r8, pcol, lchnk)
!  call outfld( 'BEAST', beast(:pcol,lchnk)*1.e-5_r8, pcol, lchnk)
!  call outfld( 'BDOWN', bdown(:pcol,lchnk)*1.e-5_r8, pcol, lchnk)
!  call outfld( 'BMAG', bmag(:pcol,lchnk)*1.e-5_r8, pcol, lchnk)
!  call outfld('WE      ', we, pcol, lchnk)

#ifdef SW_DEBUG
   if( lchnk == 25 ) then
      write(iulog,*) ' '
      write(iulog,*) '---------------------------------------'
      write(iulog,*) 'iondrift: elfld @ lchnk,i = ',lchnk,' 14'
      write(iulog,'(1p,3g15.7)') elfld(:,14)
      write(iulog,*) 'iondrift: bdown,bnorth,beast,bmag @ lchnk,i = ',lchnk,' 14'
      write(iulog,'(1p,4g15.7)') bdown(14,lchnk), bnorth(14,lchnk), beast(14,lchnk), bmag(14,lchnk)
      write(iulog,*) '---------------------------------------'
      write(iulog,*) ' '
   end if
#endif

   end subroutine iondrift

   subroutine cal_mlt( mlt, lchnk, ncol )

!-------------------------------------------------------------------------------
! Purpose: calculate the magnetic local time of WACCM geog. point
!
! Method: using the location of the geomagnetic dipole north pole,
!    the subsolar point location and the apex longitude of the 
!    geographic WACCM point the magn. local time can be calculated
!    subroutines from Roy Barnes HAO Feb. 2004
!
! Author: A. Maute Feb 2004  
!-------------------------------------------------------------------------------

   use mo_apex, only:  &
     alonm,   &  ! apex mag longitude at each geographic grid point (radians) 
     colatp,  &  ! geocentric colatitude of geomagnetic dipole north pole (deg)
     elonp	 ! East longitude of geomagnetic dipole north pole (deg)
   use time_manager, only : get_curr_calday, get_curr_date
   use cam_history,  only: outfld
   use cam_control_mod, only : magfield_fix_year
!-------------------------------------------------------------------------------
!	... dummy arguments
!-------------------------------------------------------------------------------
     integer, intent(in)   :: ncol      ! np. of atmospheric columns
     integer, intent(in)   :: lchnk     ! chunk identifier
     real(r8), intent(out) :: mlt(ncol) ! mag.local time of WACCM geo. grid point

!-------------------------------------------------------------------------------
! local arguments
!-------------------------------------------------------------------------------
     integer  :: i
     integer  :: iyear, iday, ihr, imn, imo, iday_m, tod 	! time of day [s] 
     real(r8) :: collonm, sbsllat, sbsllon, sec, ut

!-------------------------------------------------------------------------------
! get current calendar day of year & date components 
! valid at end of current timestep
!-------------------------------------------------------------------------------
     call get_curr_date (iyear,imo,iday_m,tod)  ! year, time of day [sec]
     iyear = magfield_fix_year
     iday  = get_curr_calday()                  ! day of year

!    ut   = tod/3600._r8			! UT of day [hrs]
!    ihr  = int(ut)		                ! hours
!    imn  = int((ut - ihr)*60._r8)	        ! minutes      
!    sec  = (ut - ihr - imn*60)*3600._r8        ! seconds
     ihr  = tod/3600
     imn  = mod( tod,3600 )/60
     sec  = tod - 60*(ihr*60 + imn)

!-------------------------------------------------------------------------------
!  find subsolar geographic latitude and longitude 
!-------------------------------------------------------------------------------
     call subsol( iyear, iday, ihr, imn, sec, sbsllat, sbsllon )

!-------------------------------------------------------------------------------
!  computes magnetic local time   
!-------------------------------------------------------------------------------
     do i = 1,ncol
       collonm = alonm(i,lchnk)*rtd   ! mag lons (deg)
       call magloctm( collonm, sbsllat, sbsllon, colatp, elonp, mlt(i) )
     end do

!    call outfld('MLT     ', mlt, ncol, lchnk)

     end subroutine cal_mlt

   end module exbdrift
